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We study the numerical solution of a nonlinear, partial-differential 
equation that describes charge transport in a model of a charge-coupled 
device (ccd). This model differs from previous models in that field-aiding 
of the transfer is taken into account. Although a derivation of the transport 
equation is given, the main emphasis in the paper is on the numerical 
techniques involved, and no actual numbers are presented. Actual numerical 
results based on the techniques developed here can be found in several 
recent design studies. The equation, which is parabolic, has one space 
dimension and one time dimension. Galerkin's method, with standard 
chapeau functions, is used to discretize in space. This results in a very 
stiff system of nonlinear, ordinary, differential equations. To solve these 
equations, we use a first-order backward Euler scheme coupled with 
extrapolation. A number of alternative schemes were tried and found to be 
inadequate. 

I. INTRODUCTION 

In this paper, we study the numerical solution of a nonlinear, partial- 
differential equation that describes charge transport in a model of a 
charge-coupled device (ccd). The emphasis is on the numerical tech- 
niques involved, although a derivation of the equation is given. The 
reader is referred to other papers where the solutions are used in device 
theory and design. 12 We briefly summarize the physical background 
of the equation first. 
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A knowledge of the dynamics of charge transfer in a ccd is, of 
course, central to a complete understanding of its operation. A calcu- 
lation of the motion of charge in a ccd, starting from the coupled, 
nonlinear Poisson and charge-conservation equations and taking into 
account the full geometry of the device, has so far proved impossible. 
However, Strain and Schryer 3 and, independently, Kim and Lenz- 
linger 4 developed and studied an approximate, one-dimensional model 
of charge transfer in a ccd. The original analysis considered motion 
owing only to diffusion and the mutual repulsion of the charge carriers. 
Field-aided transfer was ignored. Since these original studies, a number 
of other authors have studied the effects of field-aiding. 5-8 In Refs. 
5, 6, and 8, as in the original papers, 34 an infinite sink for the charge 
at one end of a cell is assumed. The assumption of an infinite sink 
rules out charge "bunching," which in certain situations is an im- 
portant effect (for an example of this, see Ref. 1, Fig. 8). In Ref. 7, 
the assumption of an infinite sink is not made. In this paper, we 
extend the original work 3 - 4 to include field-aiding and more realistic 
boundary conditions. Our model can describe both surface 9 ccds and 
buried-channel 10 ccds (bccds) . We do not include the effects of surface 
traps, since the main application 1 was to bccds. We feel the numerical 
scheme described here has advantages over that used in Ref. 7, where 
essentially the same model as ours was used to study surface ccds, 
with the effect of traps included. Calculations using our methods 
show that bccds, which can be fabricated with present technology, 
should be extraordinarily fast and efficient and have reasonable 
charge-carrying capabilities. Transfer times of 1.8 ns are predicted for 
a two-phase device having 10-/xm-wide electrodes. 1 Slower but similar 
results are obtained for surface devices. 

Strain and Schryer 3 solved, by the method of finite differences, a 
transport equation quite similar to the one we study here. However, 
their method of solution proved inadequate when applied to our 
equation. It is possible to obtain solutions of the transport equation 
as follows. We use Galerkin's method 11 with standard chapeau func- 
tions in space. We treat the time behavior by polynomial extrapolation 
to the limit of the results of a first-order, fully implicit (nonlinear), 
finite difference scheme. Although the equation only roughly models 
the true physical situation, an accurate knowledge of the solution as 
it varies over many orders of magnitude is necessary if it is to be of 
any use. This requirement makes the numerical solution of the equa- 
tion difficult. Many other schemes were tried, and the above method 
is the only one we found that could solve the problem. 

The equation of charge transport is derived in Section II, although 
some more complex details are given in Appendix A. The technique 
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for numerically solving the equation of charge motion is given in 
Section III, with some details in Appendix B. Questions of existence 
and accuracy are discussed in Section IV, along with the use of poly- 
nomial extrapolation. An outline of the theory of extrapolation is 
given in Appendix C. The method by which initial solutions are 
obtained is the subject of Section V. Finally, in Section VI we discuss 
several other schemes by which we tried to solve the equation of 
charge motion and which failed. 

II. DERIVATION OF THE TRANSPORT EQUATION 

We refer the reader to the literature for a discussion of the principles 
of operation of either surface ccds 9 or bccds. 10 Basically, however, both 
are devices that move packets of charge from under one electrode to 
under another electrode by suitably changing the voltage on the 
electrodes. 

As in Ref . 3, we assume that the charge can be described by a charge 
density q(x, t). Here, a; is the distance under the plates (see Fig. 1) 
and t is the time. Then, as we show in Appendix A, the component of 
the electric field along the direction of motion of the charge, which is 
due to the mutual repulsion of the charge, is 

El = - Sq x . (1) 

The elastance S is assumed to be a constant independent of x and t. 
In all that follows, we use subscripts to denote differentiation; thus, 
q x = dq(x, t)/dx, etc. Equation (1) holds for both surface and buried 
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Fig. 1 — Schematic of a ccd showing relation to device of ar-coordinate in transport 
equation. 
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channel devices, although the values of S are different in each case. 
Expressions for S are given in Appendix A in terms of the physical 
parameters of the devices. 

Let (p(x, t) be the given driving potential due to the voltages applied 
to the electrodes. For a surface ccd, <p is the electric potential at the 
oxide semiconductor interface, while, for a bccd, <p is the potential at 
the potential minimum of the buried channel. In most applications, 
we have approximated <p by the potential in the ccd in the absence of 
any mobile charge. 12,13 

The total field along the direction of motion is 

E x = - Sq x - <p x . (2) 

The current density is 9 

J(x, t) = qt iE x - Dq x , (3) 

where D is the diffusion constant and n is the mobility, which we also 
assume to be constant. If we substitute (2) into (3) and make use of 
the Einstein relation D = (kT/e)n = an, then 

J(x, t) = - M [(« + Sq)q x + q<p x ~]. (4) 

If we substitute (4) into the charge-conservation equation, 14 

qt + J x = 0, (5) 

we get the desired transport equation, 

q t = m[(« + Sq)q x + ff<pj*. (6) 

The appropriate solution of (6) satisfies an arbitrarily given initial 
distribution of charge q(x, 0) and the boundary conditions J(0, t) 
= J(L, t) = 0. The boundary conditions state that there is no charge 
flow into or out of the device at either end. L is the length of the device. 
It is convenient to write (6) in terms of dimensionless quantities, 
as in Ref. 3. Let 

t = t/(L 2 /nv ), y = x/L, w = Sq/v , $ = tp/v , = a/v„, (7) 

where v is a reference voltage. Then (6) becomes 

w T = l(w + &)w v + V&y~]y. (8) 

As it turns out, there seems to be no natural voltage unit in the 
problem (Ref. 3), so we typically pick v = 1 volt. 

Physically, the quantity of interest is the total charge present 
between any two points £> yi < y* £ 1. This suggests that, instead 
of w(y, t), we consider 

Q(V,t) = [ y w(t,T)dl (9) 

Jo 
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If we integrate eq. (8) with respect to y from to y and make use of 
the boundary condition J(0, t) = 0, we get 

Qr = (Q„ + &)Qw + <?„*,. (10) 

Since the right-hand side of (10) is just proportional to J(y, t), we 
see that Q T (1, t) = 0. From this last remark and (9), it follows that 
the correct boundary conditions on Q(y, t) are 

Q(0, r) = 0, Q(l, t) = Qt = const. (11) 

The appropriate initial condition is determined from w(y, 0) by 
setting r = in (9). The transport problem we wish to solve is, 
thus, eq. (10), subject to boundary conditions (11) and given initial 
conditions. This is a much simpler problem than attempting to solve 
(8) for the charge density. 

III. SOLUTION OF THE TRANSPORT EQUATION 

We simplify the notation slightly by setting 

Hv, r) = My, r), (12) 

and note that (10) can be written 

-0Qw -\§y (Qv) 2 - +Qv + Qr = o. (13) 

If we multiply both sides of (13) by a continuous, piece- wise differ- 
entiable function f(y) which satisfies /(0) = /(l) = 0, integrate the 
result from to 1, and integrate the terms containing second deriva- 
tives by parts, we obtain (letting /' = df/dy) 

1 ZPQv + h(Qv) 2 y(y) + L-1Q, + Qrlf(y)}dy = o. (14) 



/: 



Equation (14) is the starting point for the application of Galerkin's 
method, because any twice-differentiable function Q(y, t) that 
satisfies (14) for all continuous, piece-wise differentiable f(y) satisfying 
/(0) = /(l) = must also be a solution of (13). 

We now discretize in space by introducing a net { y\, yz, • • • , 2/at } 
on [0, 13 and a set of standard chapeau functions fj(y), 1 ^ j ^ N, 
as pictured in Fig. 2 and defined in Appendix B. In all that follows, the 
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Fig. 2 — Discretization of the space interval and the corresponding chapeau functions. 
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net {yi, • * • , 2/at } is assumed to be given and fixed. In terms of the 
basic chapeau functions, we define approximations to the solution and 
external field : 

Q(y, r) = *Z QMMv) + QrMv), (is) 

Hy,r) = Y.Ur)Si(y). (16) 

Note that Q(y, t) has been constructed to satisfy the boundary 
conditions, Q(0, t) = 0, Q(l, t) = Q T . The functions Q t (r) are yet 
to be determined, but we require that they satisfy the initial conditions 

Qy(0) = Q(yj,o). (17) 

Because of (17), Q(y, t) satisfies the correct initial conditions at the 
mesh points: Q(yj, 0) = Q{y u 0). We define ^,(r) = rf/(y jf r), so that 

UVh r) = HVii t). 

To determine the N — 2 functions Qj(t), we require that Q(y, t) 
satisfy (14) for each of the N — 2 choices of f{y), f(y) = fj(y), 
2 ^ j ^ N — I, with \f/(y, t) replaced by \f/(y, t). This yields a system 
of N — 2 first-order, nonlinear, ordinary differential equations for 
the Qj(t). This technique has a robust history and has been applied, 
not only to many problems of the same type as (13), but to other 
types of problems as well. The idea is quite simple : Let the approxi- 
mate solution be a linear combination of the functions/,-, 2 ^ j ^ N — 1 
and then make the left-hand side of (13) orthogonal to each of these 
functions. In geometrical terms, this means making the left-hand side 
of (13) orthogonal to the span of fz, ■ ■ • , /jv-i, denoted by (fa, • • • , /at-i), 
in £ 2 [0, 1] in the usual inner product: (/, g) = fo 1 f(y)g(y)dy. Then, 
crudely speaking, as more points y, are chosen, (/ 2 , • • • , /v-i) spans 
more of £ 2 [0, lj and the left-hand side of (13) must go to zero as 
N — »■» , so long as it remains orthogonal to (/2, • • •, /at-i). 

If we carry out the substitution of (15) and (16) into (14), with 
f(y) = fi(y) f 2 £ i ^ N — 1, the N — 2 equations result: 

"£ l Qi\e Pfiftdy-'t+k P f'jfkfdy + Qt f fiftfsdy] 

y=2 L •'o fc=i ^o ./o J 

+h N £ 1 Q&k f 1 f'jflcf'idy + E 1 &(r) Cfifidy 

j,k=2 Jo j=2 JO 

= - pf jj'iU'NYdy + Q T (/3 fif&du 

- t^ k f o fkfwfdy)]- (18) 
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we set Qi(m + 1) = Qi(m) + r.(ra), substitute this into (21), and 
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